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Abstract 

A new projection method based on radial basis functions (RBFs) is presented 
for discretizing the incompressible unsteady Stokes equations in irregular geome¬ 
tries. The novelty of the method comes from the application of a new technique 
for computing the Leray-Helmholtz projection of a vector field using general¬ 
ized interpolation with divergence-free and curl-free RBFs. Unlike traditional 
projection methods, this new method enables matching both tangential and 
normal components of divergence-free vector fields on the domain boundary. 
This allows incompressibility of the velocity field to be enforced without any 
time-splitting or pressure boundary conditions. Spatial derivatives are approx¬ 
imated using collocation with global RBFs so that the method only requires 
samples of the held at (possibly scattered) nodes over the domain. Numerical 
results are presented demonstrating high-order convergence in both space (be¬ 
tween 5th and 6th order) and time (up to 4th order) for some model problems 
in two dimensional irregular geometries. 
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1. Introduction 


Global Radial Basis Functions (RBFs) are increasingly used as the building 
blocks of methods for the numerical solution of partial differential equations 
(PDEs), primarily due to their ability to discretize differential operators on 
arbitrary geometries using scattered node layouts. The potential for spectral 
accuracy of global RBFs on smooth problems further adds to their appeal. 
These methods have been used to solve PDEs on planar regions [22, 21, 31, 5, 12], 
spherical regions [7, 8, 19, 36], and general surfaces [30, 17]. 


In this work, we present a novel high-order RBF collocation method for the 
numerical solution of the unsteady (or time-dependent) incompressible Stokes 
equations on some (irregular) 2D domain 12. These equations are obtained in 
the zero Reynolds number limit of the incompressible Navier-Stokes equations 
and are given by 


du 

Ik 


—-Vp -I- v/S.u -\- 
P P 

V • M = 0, 


( 1 ) 


where « = -yj is the velocity field, p is the pressure, / is some forcing term, 
p is the (constant) fluid density, and v is the coefficient of kinematic viscosity. 
We consider no-slip conditions on the boundary of 12, denoted by 5f2, 


u{x,t) = a: G (9f2, 


( 2 ) 


with the restriction that u ■ n = g ■ n = 0, where n is the respective outward 
normal unit vector on dCl. 

One of the dominant approaches to numerically solving the unsteady Stokes 
equations (and by extension, the Navier-Stokes equations) is to use so-called 
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fractional-step or projection methods, which were first introduced independently 
by Chorin [4] & Temam [33] and have advanced to more modern exact and ap¬ 
proximate projection methods based on pressure Poisson equations; for example, 
see [3, 24, 25, 38]. Broadly speaking, these methods employ operator splitting, 
and use the pressure to project an intermediate velocity field to the space of in¬ 
compressible or divergence-free velocity fields. Such methods are more efficient 
than methods that solve the saddle-point-like systems that arise from the dis¬ 
cretization of (1) (for more, see [2]). However, they typically require specialized 
grids, and the careful selection of pressure and intermediate velocity boundary 
conditions to match the actual boundary conditions on the velocity field; in¬ 
deed, these methods typically only match the normal components, and attempt 
to match tangential components through time extrapolation [3]. As a result of 
operator splitting, it is difficult (and in many cases impossible) to attain high 
orders of accuracy for the velocity field in time. In addition, errors in computing 
the pressure may propagate to the velocity field, since the pressure is often used 
to project the intermediate velocity field. 


The RBF method we develop for the unsteady Stokes equations also employs a 
projection-based approach, but avoids the issues detailed above. This is achieved 
by developing a numerical approximation to the projection operator that can 
be applied to (1) in such a way as to avoid any errors from operator splitting. 
The method relies on the following alternative form of (1) originally proposed 
by Leray [23, 9]: 



( uAu + - 

V P 



( 3 ) 


where the operator H is the so-called Leray-Helmholtz (or simply Leray) pro¬ 
jector and has the property that for a vector field w on a connected domain Cl, 
n(w) = V where V • v = 0 and v • n = 0. The existence of this operator follows 
from the Helmholtz-Hodge decomposition theorem. Eq. (3) is simply a forced 
diffusion-type equation for u and is obtained by applying H to the first equation 
in (1) and using the fact that H (Vp) = 0. Using the orthogonal complement 
of n, which we denote as H^, the following auxiliary equation to (3) for the 
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gradient of the pressure can be obtained: 

-Vp = U^ (i^Au+-f] . (4) 

P V P J 

This allows the pressure to be recovered as an instantaneous functional of the 
velocity field. Our goal is to construct a discrete approximation to 11 and solve 
(3) using collocation and the method-of-lines. We will use this to recover the 
pressure as in (4). 

To discretize the Leray projector, we adopt the recent method from [18] for 
computing the Helmholtz-Hodge decomposition of a vector held based on gen¬ 
eralized RBF interpolation with matrix-valued kernels. This method provides 
a well-posed way to construct a discrete Leray projector over a set of scattered 
nodes X in some domain Cl, with the exact normal boundary conditions enforced 
at points along the domain boundary. It can then be applied to any vector held 
sampled at X to recover an analytically divergence-free held everywhere on in¬ 
terior and boundary of Cl. We modify the technique of [18] here so that 11 also 
incorporates tangential boundary conditions, unlike both the traditional Leray 
projector and the implicit Leray projection performed by methods that use the 
pressure Poisson equation. An approximation to 11^ can also be obtained triv¬ 
ially from this technique, and we use that to recover the pressure p, but without 
solving a differential equation for p . We combine the discrete Leray projection 
with a global RBF collocation method for approximating the Laplacian in (3) 
to get a semi-discrete approximation, which we then integrate in time with a 
high-order BDF scheme. The method avoids any errors from operator splitting, 
does not require the use of specialized grids or meshes, and can be used with 
scattered nodes. We demonstrate that this new methodology allows for high 
order accuracy in space and arbitrary orders of convergence in time on irregular 
domains. 

Similar approaches to the one we propose have been adopted with divergence- 
free RBFs [35], divergence-free wavelets [20], and divergence-free polynomi¬ 
als in the context of hnite elements [37]. The first approach uses a simi- 
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lar matrix-valued kernel, but enlarges it to account for the pressure coupling, 
which increases the computational cost. It also is presently restricted to time- 
independent problems. The wavelet approach is similar to ours in that it em¬ 
ploys the (standard) Leray projection, but it is restricted to far more regular 
domains and nodes (rectangular boxes with tensor product nodes). The finite 
element approach differs from ours in that a weak formulation with divergence- 
free test functions is employed to remove the pressure. This method also requires 
a mesh, and boundary conditions are only enforced weakly. 

The remainder of the paper is organized as follows. In the next section, RBF col¬ 
location for approximating the Laplacian is briefly reviewed. In Section 3, gener¬ 
alized interpolation of vector-valued functions with matrix-valued RBFs is pre¬ 
sented and applied to the construction of a discrete Leray projector. Equipped 
with spatial discretizations of the Laplacian and the Leray projector, we then 
discuss the full space-time discretization of the incompressible unsteady Stokes 
equations in Section 4. Eigenvalue stability of the method is investigated in Sec¬ 
tion 5. Numerical results demonstrating the spatial and temporal accuracy of 
our method are given in Section 6 using two problems involving time-dependent 
boundary conditions: rotational flow on a rotating disk, and a flow on a rect¬ 
angular domain containing a spinning disk in its interior. This is followed by a 
summary of the method and a discussion of future enhancements in Section 7. 


2. RBF Collocation method for the Laplacian 

Let fl C and >]Rbea kernel with the property <j){x,y) := 

4>{\\x — y\\) for x,y & fl, where || • || is the standard Euclidean norm in We 
refer to kernels with this property as radial kernels or radial functions. Given 
a set of nodes X = C f2 and a continuous target function / : 12 —> M 

sampled at the nodes in X, the standard RBF interpolant to the data has the 
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form 


N 

®/(®) = X^Cfe(/)(||a: - XkW). 
k=l 

The expansion coefficients {ck}^^i are determined by enforcing 
This can be expressed as the following linear system: 


(^(ri,2) . 

4’{r2,i) (t>{r2,2) ■ 

■ 4’{r2,N) 


Cl 

C2 


II 

f2 

HrN,i) <P{rN,2) ■ 



Cpf_ 


Jn_ 


Cf fx 


(5) 

f\x- 

( 6 ) 


where rij = \ \xi — Xj\\. li (f) is, for example, a positive-definite radial kernel on 
M'^, and all nodes in X are distinct, then the matrix Ax above is guaranteed 
to be positive definite, hence (5) is well-posed. Examples of various choices for 
<j), including relaxed conditions to guarantee the well-posedness of (5) can be 
found in [5, Ch. 3-12]). The convergence of these RBF interpolants to the target 
function as the number of samples increases is, in general, determined by the rate 
of decay of the Fourier transform of (p- If the rate of decay is algebraic, then one 
can expect algebraic convergence for a sufficiently smooth function. If instead 
the Fourier transform decays exponentially, then one can expect exponential (or 
spectral) convergence for a large class of infinitely smooth functions. 


The RBF interpolant (5) can be exploited to construct discrete approximations 
to linear differential operators (differentiation matrices) in the same fashion as 
standard Chebyshev or Fourier collocation methods [12, Ch. 3-4]. However, the 
RBF method is not restricted to special node sets. Let £ be a linear differential 
operator of interest and suppose we wish to approximate Cf using the RBF 
interpolant (5) of / sampled at X. Applying C to (5) and evaluating at the 
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points in X gives the following matrix equation: 


{Lsf)i 

{LSf)2 


jC<j){ri^i) jC(p{ri^2) ■ 

C(f){r2,i) L(t>(r2,2) • 

• C(l>{ri^N) 

■ L4>{r2,N) 


Cl 

C2 

iCSf)N_ 


£(^(rjv,i) L(j){rN,2) • 

. C(f>{rN,N) 


cn 



where {Csf)j = and Ctpirj^k) = £<^(||a; - Prom (6) we 

know that c/ = A'j}fx, so that the right hand side of (7) can be simplified to 
Aj^A^^fx- The matrix 


Lx ■— Li^A^, 


( 8 ) 


gives a discrete approximation to the operator C over the node set X and is 
typically called a differentiation matrix for C. If s/ gives a good approximation 
to / then we expect Lxfx to give a good approximation to Cf at X. In this 
study, we are interested in £ = A (the Laplacian) which makes the entries ^4^ 
simple to determine since the Laplacian in of a radial kernel <j) is Xcj) = 
where the derivatives are taken with respect to the univariate 
radial variable r. Boundary conditions can be implemented in a similar fashion 
to traditional spectral collocation methods [34]. For Dirichlet conditions (which 
arise for the assumed no-slip boundary conditions in this study), it makes sense 
to partition the nodes in the set X before computing the differentiation matrix 
(8) so that the interior nodes, Xj, of Q, appear first, followed by the boundary 
nodes Xb- If Nj and Nb denote the respective numbers of nodes in Xj and 
Xb, then we can implement these boundary conditions by simply removing the 
last Nb rows of Lx- We denote this new matrix as Lx and partition it further 
as 


Lx = [Li Lb] , (9) 

where Li is an Ni x JV/ matrix and Lb is an Ni x Nb matrix. The matrix 
Li operates on samples of a function at the interior nodes and Lb operates on 
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samples at the boundary nodes. For a more detailed discussion of differentiation 
matrices and boundary conditions for RBFs, see [5, Ch. 42]. 


3. Approximating the Leray Projector 

We now focus on our approach for discretizing the Leray projector 11. Recall 
that if / is a vector field on a bounded domain Cl with Lipschitz boundary dCl 
then / can be uniquely decomposed as / = n/ + 11^/, with 11/ divergence- 
free and tangential to dCl, and 11^/ = Vg for some scalar potential q [9]. 
In the popular class of projection methods, q is either an approximation to 
the fluid pressure, or is closely related to it (depending on the precise form of 
the projection method used) [3]. While traditional (fractional-step) projection 
methods implicitly apply the Leray projection by calculating q and subtracting 
off Vg from /, we seek to approximate 11 explicitly using only samples of the 
field / at scattered nodes. To accomplish this, we use the vector decomposition 
method outlined in [18] based on generalized interpolation with divergence-free 
and curl-free kernels. 

3.1. Divergence-free and Curl-free Kernels 

Divergence-free and curl-free matrix-valued kernels have been well-studied, both 
on planar domains and on the sphere [27, 26, 28, 14, 16]. On these kernels 
are given by 

^div = (-A/ -I- VV^) ((> and ^curi = -VV^((>, (10) 

respectively, where (j) is a, standard scalar-valued RBF. Each kernel is matrix¬ 
valued (a map from x x M'^), and the columns of ^div and ^curi are 

divergence-free and curl-free, respectively [18]. Letting r = [[a:— a:jH, x = 
and '0 = yX0r), these kernels are given explicitly in d = 2 dimension for inputs 



X = (x,y) and Xj = {xj^yj) as 


^div ) 


^curl (^1 ) 


-X{r) - 'ipir){y - yj)'^ 
^p{r){x - Xj){y - yj) 

-X{r) - il}{r){x - XjY 
-xp{r){x - Xj){y - yj) 


i!{r){x - Xj){y - yj) 
-X{r) - ip{r){x - Xjf 

-^{r){x - Xj){y - yj) 
-X{r) - 'tp{r){y - Vj)^ 


and (11) 


( 12 ) 


The second argument in each term acts as a shift of the kernel in the same way 
that it does in the scalar RBF setting. 


Another related kernel, which we refer to as the “full” kernel, is given by 


^ = ^div + ^curi = 


where / is the dxd identity. This diagonal kernel splits naturally into divergence- 
free and curl-free parts. In what follows we will let ^div '■ L 2 (IR'^) —> L 2 (IR'^) 
denote the projector to the subspace of divergence-free fields in L 2 (ffi'^). It can 
be shown that = ^div, where fdiv is applied to the columns of and 

that $ = ^div + ^curi gives the L 2 (K'^) Helmholtz-Hodge decomposition of $ 
(interpreted column-wise) [18]. Even though this decomposition is independent 
of Cl and hence does not yield the correct boundary conditions with respect to 
n, we will still make use of it in deriving the kernel-based approximation to the 
Leray projection. 

One can construct a generalized interpolant to a vector held by linearly com¬ 
bining shifts of columns of one or more of the basis functions mentioned above. 
The form of the approximation depends on the requirements one wants to im¬ 
pose. For example, if one wants to interpolate a vector held at a point Xj, a 
term of the form ^{■,Xj)cj should appear in the approximation, where Cj is a 
d X 1 vector of coefficients. If one wants to require that the divergence-free part 
take a certain value at a point yj, a term of the form ^div{xyj)^j should be 
included. For more details, see [18]. 
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3.2. Kernel-based Leray Projector 


Let Xj = {xk}kLi be the set of nodes on the interior of the domain Vl, and 
= {yj}f=i be the set of nodes on the domain boundary dfl. We will ap¬ 
proximate n/ by constructing a generalized interpolant that matches / on the 
interior nodes Xj and matches the correct divergence-free boundary conditions 
at the boundary nodes Xb. We note that even though the Leray projector is 
completely determined by specifying the normal component of the divergence- 
free term, we can also impose tangential boundary conditions if they are avail¬ 
able to usd The generalized interpolant we use is of the form 

Ni Nb 

{■,yj) dj, (13) 

fc=l j=l 

where Cj and dj are d x 1 vectors of unknown coefficients. To solve for the 
coefficients, we impose the conditions 


^divSfiy^) = siy^), 


£ =1,.. .,Ni, 
m = 1,.. .,Nb, 


(14) 

(15) 


where g represents the given boundary conditions for 11/. This leads to a 
symmetric linear system of the form: 


A 

B 


c 


fl 

B^ 

C 


d 


gB 


(16) 


where c and d contain the coefficients in (13), f ^ and gs represents the data 
coming from f on the interior and g on the boundary nodes, respectively, and 
the matrices A, B and C arise from applying the conditions in (14) -(15) to the 
basis functions in (13). Specifically, A, B and C are given in d x d blocks as 


^This is in contrast to traditional approaches that recover the Leray projection, which 
typically are only able to enforce normal boundary conditions and show that the tangential 
component is within some truncation error [3]. 
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follows: 


^{xi,Xl) ■■■ $(£C1,£CJV/) 


A = 


B = 


C = 


^(xNj,Xl) ■ 
^div{xi,y^) 

^div{xN,,yi) 

^div{yi,yi) 


^{xNr,XNi) 

^divixi,yN 


^div (xNi I y Nb ) 

^div{yi,yNB) 


^ div (yjvj>yi) ••• ^ div [yNi^yNs) 


(17) 


(18) 


(19) 


The {dNj + dNs) x {dNj + dNs) matrix in (16) is symmetric and positive- 
definite [18]. Thus this system is always invertible, even for scattered nodes X. 
We remark that in this work, d = 2- however, this approach naturally extends to 
the case of d = 3 as well. Error estimates for generalized interpolants like (15) 
generally depend on the decay rate of the Fourier transform of the underlying 
radial kernel cj) as in the scalar RBF case [18]. 


Note that the definition of A, B and C above coincides with the standard form 
of the scalar RBF matrix given in (6). Written in this way, the rows and columns 
of each are interlaced in an alternating x coordinate, y coordinate pattern. From 
an implementation perspective, it is more convenient to organize the matrices so 
that the data is given in coordinate blocks. Since the evaluation of each kernel 
is a 2 X 2 matrix, we organize A, B and C into four blocks each, where each 
upper left block corresponds to the evaluations of the upper left kernel function 
(cf. (11)), each upper right block corresponds to the evaluations of the upper 
right kernel function, etc. Denoting the re-ordered matrices as A, B, and C, 
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(16) becomes 


A B 


c 

B'^ C 


-1 

_1 


ff 

fi 

9% 

9% 


( 20 ) 


where the data vectors have now been rearranged so that ff and ff represent 
the X and y coordinates of f evaluated at Xj, and gg and Qg represent the x 
and y coordinates of g evaluated at Xg. Note that due to the diagonal structure 
of the full kernel this has the effect of transforming the 2Nj x 2Nj matrix A 
into a block diagonal matrix A, where the (i,i)^^ coordinate of each identical 
Nj X Ni diagonal block in A is given by —X<f){\\xi — Xj\\), i,j = 

Since in practice the bulk of the entries are taken up by A, storing the matrix 
in this form reduces storage costs. Further, the block diagonal structure makes 
it possible to solve (20) using a more efficient Schur complement method than 
if the upper left portion was left interlaced. 


Once (20) is solved, an approximation to the Leray projection of / on Xi is 
obtained by simply evaluating the divergence-free part of sy, which is given by 

Ni Nb 

^div^f ^ ^ ^div (* 5 ^ ^ ^div (* 5 U j ) • 


fc=l j = l 

To this end, we need to evaluate shifts of the divergence-free kernel ^div on the 
interior nodes. Let A^iv be the 2Ni x 2Ni matrix obtained as in (17), but with 
the full kernel replaced by ^div, and Adiv denote the reordered Adiv according to 
the permutation of A described above. With this, the discrete Leray projector 
P which approximates 11 on X^ is then given by 

r ^ n -1 


Px = Adiv B 


A B 
C 


= P, 


( 21 ) 


where we have partitioned Px into the 2Ni x 2JV/ matrix Pj and the 2Nd x 2Nb 
matrix Pb, which operate on interior and boundary samples, respectively. 
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3.3. Recovering the Pressure 


In this section we illustrate how to extract from an approximation to the 
scalar potential q, where = Vg. First, we consider the curl-free part of Sf, 
which is given by Vcuri^f = ^div^f consists of shifts of the curl-free kernel 
from (10). Expressing the kernel in terms of derivatives of (j) as in (10), one can 
factor out a gradient to obtain 

N, Ni 

^curl^f — ^ ^ ^curl (*; ^fc) ^ ^ ^ ^ 1^(11* ^/^ll) ^k t (^^) 

k=l k=l 

V 

where the scalar potential ip in (22) approximates q (up to an additive constant). 


A discrete operator for approximating q at the interior nodes Xj can be con¬ 
structed in a similar fashion to how the discrete Leray projector was con¬ 
structed in (21). To this end, we define matrices and Ay with entries 
{Ax)i,j = {xi - Xj)x{\\xi - XjW) and {Ay)ij = {tji - yj)x{\\xi - Xj\\), for 
z, j = 1,..., Nj. The discrete operator Qx that evaluates the gradient potential 
of n^/ is then given by 



A B 
B'^ C 


Qi Qb ■ 


(23) 


4. Method of Lines Formulation 

We use a method of lines formulation to solve the incompressible unsteady 
Stokes equations (3) with boundary conditions (2) using the discrete versions 
of the Laplacian (9) and Leray projector (21) presented in the previous two 
sections, respectively. As in those sections, we consider a set of nodes X over 
the domain fl and partition them into a set of Nj nodes. A/, on the interior of 
Q, and a set of Nb nodes, Xb, on the boundary dfl. 

Let ui, vj, /f, and ff denote vectors containing samples at Xj of the u and 
V components of the velocity vector u, and the x and y components of the 
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forcing function /, respectively, in (3). Similarly, let and denote vectors 
containing samples of the x and y components of the boundary conditions g 
in (2) at Xb- Using the discrete approximations to the Laplacian and Leray 
projector, we discretize (3) in space at the nodes Xi and arrive at the following 
system of ordinary differential equations: 


d 

dt 


Ui 

Vi 


( 




_ _ 




_ _ 



Lj 0 


Ui 

+ v 

Lsgl 

Pi 

ly 






0 Lj 


Vi 


LBgl_ 


1 ^ V ^ 





\ Li 


\ 



_ _ 



_ _ 

1 

-b - 

ff 


+ Pby 

9% 

p 

Jl 


dt 

/b_ 







) 

(24) 


The fact that time derivatives of all components of the field on the boundary 
are included in the last term deserves some discussion. The unsteady Stokes 
equations (1) can be written in the form of the Helmholtz-Hodge decomposition 
theorem: 

« 1 du 1 , ^ 

+-/=—+-Vp, (25) 

p p 

W ^ Vq 

where V • v = 0 and n • v = 0. The latter condition comes from the boundary 
condition (2), with assumed restriction n • g = 0. The Leray projection applied 
to w recovers v, implying that the boundary conditions on v (given by time 
derivatives of g) can be substituted for the boundary conditions on w. As 
mentioned in the previous section, the condition n-v = 0 is enough to completely 
determine the continuous Leray projector. However, the boundary condition (2) 
also tells us what the tangential component, t-v, is on the boundary and thus all 
components are included in the RBF formulation of the discrete Leray projector 
( 21 ). 

Provided the coupled system of ODEs (24) is stable (see the next Section), it 
can be integrated in time using a suitably chosen time-integrator. In this study, 
we use backward differentiation formula (BDF) methods for this task as the 
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system is stiff. When the forcing function is available in closed form and is 
independent or linear in u, it can be incorporated easily into BDF schemes. In 
the case that it is non-linear or not available in closed form, semi-implicit BDF 
schemes would be more appropriate [1]. For BDFl (backward Euler), the fully 
discrete system for approximating (3) is given by 


I 0 
0 I 


AtPr 


— h'AtPf 


Li 

0 


0 

Li 


n+1 


..'H. + l 


( 

LBi9%)^+^ 

1 

(J.)n+1 

1 +Pb 


- (5|)" 


-b - 


\ 

LBi9l)^+\ 

P 

_(/!)"+\ 

J 


- (5^)". 


(26) 


where superscript n denotes the value of the variable at time and superscript 
n + 1 denotes the value at time -I- At. Extensions of the scheme to higher 
order BDF schemes should be straightforward. Note that, for all choices of time 
integration schemes, the approximation to the time derivative on the left hand 
side of (24) should also be applied to the time derivative in the last term on the 
right hand side of this equation for consistency. 


While we do not explicitly solve for the gradient of the pressure or the pressure 
at any time-step, these quantities can be easily recovered. An approximation to 
the gradient of the pressure can be recovered after computing and 
using a discrete approximation to (3). Specihcally, for the scheme (26) we have 


Vp 


n+11 


\Xi 




L,<+i+Lb(5|)”+i 

Liv^+^ + LB{gir+^ 


^fpn+l 


I 

iff) 


V'ln+l 
I 


J>_ 

At 


- V- 


X 


where p = pv being the coefficient of dynamic viscosity; this is just the discrete 
approximation to (25). The pressure, on the other hand, can be approximated 
(up to a constant) using the discrete operator Qi defined in (23) applied to the 
discrete approximation to w in (25). Specifically, at time-level n + 1 we can 


15 



approximate the pressure as 


P 


n+1 



\ [liv^+^ + LB{g%r+\ 


+ 


(y.)„+l 

(/y)n+l 



( 51 )”+' 

( 51 ,)”+' 



+ 


( 27 ) 


It should be clear from the above equation that the pressure is recovered purely 
as a function of the velocity at any time level. Our method does not employ 
any time-splitting. 


The matrices L/, Lb, Pi, and Pb in (24) need only be computed once as a 
preprocessing step. We solve (8) to compute Li and Lb, and (21) to compute Pj 
and Pb ■ If the pressure is also to be computed, then Qj needs to be constructed 
by solving (23), which involves the same linear coefficient matrix as (21). For 
all of these computations we use Gaussian elimination. Since these matrices 
are dense this requires 0{N^) operations up front. To solve the implicit linear 
system in the full discrete equations, such as (26), we first compute the LU 
decomposition of the system (i.e., the matrix on the left hand side of (26)) 
in a preprocessing step, which again requires 0{N^) operations. We then use 
this decomposition to solve for and u^+' at subsequent times. Each time- 
step (following any required steps for bootstrapping the time integrator) then 
requires 0{N'^) operations. This relatively high computational cost is offset 
by the fact that the method does not suffer from any splitting errors and, as 
demonstrated in the numerical results section, generates high-order accurate 
solutions in both space and time for moderate size N. In Section 7 we make 
some remarks on how the computational cost can potentially be reduced. 


5. Stability of the Numerical Method 

In this section, we investigate the time stability of our method by exploring 
the spectra of the spatial differential operator L/ and its projection PjLi in 
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(24). A necessary condition for stability is that the eigenvalues of this latter 
operator must lie within the stability domain of the time integrator used. This 
usually requires, at the very least, that all eigenvalues have non-positive real 
parts. Unfortunately eigenvalue stability is not guaranteed for general node 
layouts, but in practice we have observed that the method is stable and can be 
quite robust even with respect to node configurations that are quite scattered. 
Below we present the results of three experiments that illustrate the results 
of a much larger investigation we carried out on eigenvalue stability. One of 
these does show that stability can be an issue when the boundary is not well 
resolved. This indicates that more work needs to be done to determine sufficient 
conditions for eigenvalue stability. 

The following examples involve nodes in a domain Q, obtained by deleting a disk 
from the interior of the rectangle [0,1.5] x [0,1] (see Figure 1). This domain is 
also used one of the convergence studies in the next section. In each example, 
we used positive-definite Matern radial kernel [5, §4.4] 

<j){r) = + 15(£r)"‘ -I- 105(er)® -|- 420(£r)^ -I- 945(£r) -I- 945)e“^’', (28) 

with fixed shape parameter £ = 10 to generate both the discrete Laplacian and 
Leray projector. Similar results were seen for other positive definite kernels. All 
node sets used in the experiments can be obtained from [15]. 

The hrst example has a node conhguration that leads to unstable eigenvalues. 
We obtained the interior nodes by generating Halton points [5, §A.l], on the 
square [0,1.5] x [0,1.5] and selecting those within Cl. For the boundary nodes, we 
used 492 uniformly spaced nodes, then removed 5 from the upper-left corner; see 
the left plot in Figure 1 (a). The spectrum of Lj in this case has one positive 
real eigenvalue and PjLj has two positive real eigenvalues as shown in right 
plot of Figure 1 (a). However, as can be seen in Figure 1 (b), reinserting the 5 
points at the corner resulted in eigenvalues with all negative real parts. As a 
further test of the robustness of the eigenvalues, we began with the Halton set 
in Figure 1 (b) and removed 200 randomly selected points along the boundary. 
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(b) Halton nodes with uniformly spaced boundary nodes 



(c) Halton nodes with 200 randomly selected nodes removed from boundary 


Figure 1: The right column shows the eigenvalues of the associated spatial operators L/ and 
PjLj for the corresponding nodes in the left column. In (a) PjLj has two positive real 
eigenvalues, the largest being equal to 58.12, which is difficult to see in the figure due to the 
scale. In (b) and (c) all eigenvalues are in the left half plane. 
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This experiment was run several times and, suprisingly, we never encountered 
eigenvalues with positive real part. An example from these experiments can be 
seen in Figure 1 (c). 




Figure 2; Quasi-uniform nodes with smoothly increasing density near the boundary (left) and 
the associated eigenvalues of the associated spatial operators Lj and PjLj (right). 

The last node set considered in this section was generated in MATLAB using 
the distmesh package [29], and contains Nj = 1790 nodes on the interior and 
Nb = 492 nodes on the boundary. These nodes are quasi-uniform and increase 
in density near the boundary (see Figure 2). All eigenvalues of Lj and PjLj 
in this case lie in the left half plane, and are more tightly clustered around the 
negative real axis (see Figure 2) than the previous examples. Nodes such as 
these are utilized in the next section to test the convergence of our method. 


6. Results 

We present the results of two tests conducted with our method. In the first 
test, we compare it to an analytically known solution on a spinning disk. In 
the second test, a more irregular domain is considered and we give results from 
a refinement study on our method. In all tests we used the Matern kernel 
given in equation (28) to generate both the discrete Laplacian (9) and Leray 
projector (21). This kernel is positive definite in and its two-dimensional 
Fourier transform decays algebraically at a rate of 13. Based on the results 
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in [18], we should expect the error in the generalized interpolant (13) (and 
hence the discrete Leray projector (21)) to decay algebraically like 0(h^®) for 
sufficiently smooth vector fields when measured in the norm, where hx is a 
measure of the spacing the the nodes. Error estimates for Laplacian with this 
Matern kernel kernel are one order less, i.e., like 0{hj^^). However, for the test 
problems considered, we always observed higher rates than 4.5. It is possible 
to use two different kernels, one for the Leray projector and a smoother one for 
Laplace operator, so that the theoretical rates of convergence match. However, 
in our experiments, this made very little difference in terms of the total error 
achievable. We therefore only include results based on a single kernel. 

For both tests, we used the distmesh package [29] to generate the node sets for 
the experiments. We used the option to create a variable density node distribu¬ 
tion, with higher densities of nodes near the domain boundaries. Clustering of 
nodes at the boundaries is known to prevent Runge-type boundary effects for 
RBF methods [10]. One could also use other algorithms to generate such node 
sets [11]. 

As mentioned previously, we used BDF methods for integrating the semi-discrete 
system in time. In the first test problem, the temporal convergence of our 
method is compared for BDF schemes of orders of accuracy two to four (i.e. 
BDF2-BDF4). Lower order BDF schemes were used to generate the starting 
values for the higher order schemes in a consistent manner. The BDF4 scheme 
was used for all spatial convergence studies. Errors were computed using the 
standard £2 and foo norms. 

The node sets used in all tests, and an implementation of the method using 
BDF2 can be downloaded from [15]. 

6.1. Rotating disk with forcing: analytical solution 

This test considers the incompressible unsteady Stokes equations with forcing 
on a spinning unit disk. The forcing term f was generated so that the analytical 
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solution to the problem for all time is given as 


-7r2/sin (|(a;2 + sin(7rt) 

ii(x, y, t) = 

7ra:sin + y^)) sin(7rt) 

P{x,y,t) = sm{x - y + t). 


(29) 

(30) 


The forcing is simply f = ^ — + Vp. Boundary conditions for u were 

obtained from the analytical solution, which is purely tangential to the boundary 
and non-constant in time. For all tests the coefficient of kinematic viscosity was 
set to 1 / = I and the discrete equations were solved to a hnal time t = 1.5. Figure 
3 shows the numerical solution at t = 1.5 together with the discretization nodes. 
The streamlines shown were computed using the generalized interpolant (13) in 
an analogous technique to extracting the pressure in (22); see [18] for details. 



Figure 3: Numerical solution to the rotating disk with forcing example at t = 1.5. The solid 
lines are stream lines of the solution and arrows indicate the velocity field at a few locations 
in the domain. Solid circles mark the N — 1073 nodes used for computing the solution. 


Figure 4 shows the results from a temporal convergence study for the velocity, 
using a node set of size N = 3236. In these simulations the shape parameters 
for both the standard RBF interpolant and the generalized interpolant were set 
to e = 11. Errors in the velocity field were computed against the analytical 
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Timo Onni/ornarma nf Kliimarinal Qoliifrinn 


Time Convergence of Numerical Solution 



Time-step (At) 



Time-step (At) 


Figure 4: Comparison of the temporal convergence for the rotating disk problem using BDF2, 
BDF3 and BDF4 as the time integrators and N — 3236 nodes. The figure shows the relative 
£2 (left) and ioo (right) errors at time t — 1.5 for the velocity. The dashed lines have slopes 
corresponding to second-order {p — 2) and fourth-order [p — 4) convergence. The leveling off 
of the error is due to spatial errors dominating the temporal errors. 


solution (29). Figure 4 shows that our method demonstrates the correct order 
of convergence in time for the velocity field based on the choice of time integra¬ 
tor: second order for BDF2, third order for BDF3 and fourth order for BDF4. 
Convergence in all cases eventually levels off as the spatial errors dominate. 

Interestingly, the error for the discrete pressure obtained from our method seems 
to be unaffected by the temporal errors in the velocity. For the time-step ranges 
tested on this problem, the error in the discrete pressure appears to be de¬ 
termined purely by the spatial discretization, regardless of the time integrator 
used. 

Figure 5 (a) and (b) show the results of a spatial convergence study for both 
the velocity and pressure, respectively. The shape parameter was again fixed 
at e = 11, and the time-step in these simulations was fixed at At = 1.5/8192 
so that spatial errors dominate. Numerical solutions were computed using node 
sets of size N = 163, 369, 669, 1073, 2025 and 3236. Relative errors were 
computed against the analytical solutions (29) and (30), and are displayed in 
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Figure 5. The errors are plotted against '/N since its reciprocal gives a rough 
measure of the spacing between the nodes. Lines of best fit to the data are also 
included in the figure with the slopes reported in the legends. Based on these 
slopes, we see that part (a) of this figure shows that the velocity appears to 
converge in space at a rate between 5.5 and 6 in both the i^c and £2 norms. 
Part (b) of this figure indicates that the pressure converges at about this same 
rate in the £2 norm, but at a slower rate of 5 in the £qo norm. 

Q Spatial Convergence of Numerical Solution 

10'’ 

10'^ 

10-'' 

10" 


15 20 25 30 35 40 45 50 55 

\fN 

(a) Velocity 

Figure 5: Spatial convergence for the rotating disk problem at time t — 1.5. Figure (a) shows 
the relative £2 and ^oo errors in the velocity and (b) similarly shows the errors in the pressure. 
The dashed/dotted lines are for where C is a constant, determined by the line of 

best fit to the corresponding data. 


Spatial Convergence of Numerical Solution 



(b) Pressure 



6.2. Rotating disk in a stationary box: refinement study 

This test considers the incompressible unsteady Stokes equations in a rectangu¬ 
lar domain with a circular hole. Letting flu = [0,1.5] x [0,1] denote the rectangle 
and flc = {{x, y) : (x — 0.3)^ + {y — 0.5)^ < 0.04} the circular hole, the domain 
is given as 12 = flu\flc. To generate motion, the following time-dependent. 
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(b) Pressure contours 


Figure 6: Numerical solution to the spinning disk-in-a-box example at £ = 1.5. (a) The solid 
lines are stream lines of the numerical solution and arrows indicate the velocity field at a few 
locations in the domain, (b) The solid and dashed lines are contours of the pressure at 21 
equally spaced values over [—20,20]. The horizontal line through the center marks the zero 
contour, dashed lines indicate negative values, and colors indicating the magnitude of the 
contour arranged from darker (smaller) to lighter (larger). In both plots solid circles mark 
the N — 3077 nodes used for computing the solution. 

purely tangential boundary conditions were applied to dQc’- 

5{y — 0.5) sin (vrt) 

—5(x — 0.3) sin (nt) 


{x,y)ednc, (31) 


u{x,y,t) 
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which corresponds to a spinning disk with speed 1 at t = 1.5. At the outer 
boundary Cljj the velocity was set to zero for all time. The presence of corners 
in the outer boundary presents numerical challenges as this may add weak sin¬ 
gularities (that are damped in time) to the solution and lead to larger errors 
near these regions. 

We are not aware of an analytical solution to the incompressible unsteady Stokes 
equations in this domain with these boundary conditions. Therefore, to test the 
spatial and temporal convergence of our method we performed a refinement 
study where lower resolution numerical solutions are compared to a high res¬ 
olution solution. All simulations were run to a final time of t = 1.5 units. 
The numerical solutions for the velocity and pressure from one of these sim¬ 
ulations are displayed in Figure 6 (a) and (b), respectively, together with the 
discretization nodes that were used. Streamlines were computed as described in 
the previous section. Only results for the BDF4 time integrator are presented 
as similar results were observed for the lower order BDF methods. 



Figure 7: Temporal convergence results for the velocity field at t = 1.5 for the spinning disk- 
in-a-box example using BDF4 as the time integrator and N = 10825 nodes. The dashed line 
has a slope corresponding to fourth-order [p — 4) convergence. The leveling off of the error is 
due to spatial errors dominating the temporal errors. 


For the temporal convergence test, a set oi N = 10825 boundary-clustered 
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nodes was used, which is large enough so that temporal errors dominate in the 
range of time-steps considered, and the shape parameter for both the standard 
RBF interpolant and the generalized interpolant was set to e = 32. Errors 
were computed against the numerical solution obtained with a time-step of 
At = 1.5/2048. Figure 7 displays the results for the velocity field from the 
test and shows that the method converges to the high resolution solution at a 
slightly faster rate than fourth order. 

For the spatial convergence test, the time step was set to At = 1.5/2048 for 
all simulations. The numerical solution computed with N = 10825 nodes and 
£ = 32 was used to compute errors in the lower resolution solutions. Two types 
of convergence studies were performed. For the hrst, the shape parameter was 
set to £ = 32 for both the standard RBF interpolant and the generalized RBF 
interpolant in all simulations. With this strategy, the condition numbers of the 
interpolation matrices in (6) and (16) are allowed to steadily increase with N 
until maximum condition numbers of approximately 10^^ are reached at the hnal 
node set. For the second convergence study, we kept the condition numbers of 
the interpolation matrices roughly fixed at around 10^^ by allowing £ to increase 
with N. This strategy is motivated by the fact that while estimated convergence 
rates for RBF-based methods only hold for the fixed e case, the highest accuracy 
for a given N can typically only be obtained using the largest condition number 
that still results in numerically stable discrete operators (see, for example, [22]). 
Thus, convergence studies that use this approach are often conducted in the 
literature. 

Figure 8 displays the results from these two convergence tests. Part (a) and (c) 
of the figure show the fixed e results for the velocity and pressure, respectively, 
together with the lines of best fit to the errors. Based on these lines, we see 
the velocity errors converge at a rate between 5 and 5.5 in the £2 norm and at 
approximate rate of 6 in the foo norm. The errors in the pressure from part 
(c) show a similar convergence rate in the £2 norm, but the rate for the foo 
norm has dropped by 1. These results are not too different from the spatial 
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(a) Velocity, fixed e 



(c) Pressure, fixed e 



(b) Velocity, variable e 



(d) Pressure, variable e 


Figure 8: Spatial convergence for the spinning disk-in-a-box problem at time t — 1.5. Figures 
(a) and (c) show the £2 and £00 errors in the respective velocity and pressure using the fixed 
e approach, while (b) and (d) similarly shows the errors for the variable e approach. The 
dashed/dotted lines are for where C is a constant, determined by the line of best 

fit to the corresponding data. 


convergence study from the previous test. The velocity and pressure appear to 
be very smooth (see Figure 6), which explains why convergence is not adversely 
affected by the presence of corners. Parts (b) and (d) of Figure 8 show the 
convergence rates for the velocity and pressure, respectively, with the variable 
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e strategy. The figure shows that the convergence rates have slowed over the 
fixed £ strategy. However, for a given N (except the last N for the pressure) 
the errors are lower than the hxed e case. 

We note that these results are not entirely unexpected. Indeed, it has been 
observed many times in the RBF literature that fixing the condition number of 
the interpolation matrices (sometimes called “stationary interpolation”) results 
in a scheme that initially converges at a slower rate, but then convergence even¬ 
tually stagnates [5, Ch. 15-17]. Fixing e with increasing N leads to theoretical 
convergence, but the interpolation matrices then become more ill-conditioned. 
This however can be overcome in certain cases by using a “flat” algorithm such 
as RBF-QR [13, 6]. These algorithms have yet to be extended to the generalized 
interpolation technique presented in this article. 


7. Summary 

In this work, a new method based on the Leray projection was presented to 
solve the unsteady Stokes equations on irregular domains. The accuracy of the 
new method was demonstrated on problems with no-slip boundary conditions 
on a domain with a smooth boundary and on a domain with corners. Our 
conclusions are as follows: 

• Unlike current Leray-type projection methods, the method allows one to 
build both tangential and normal boundary conditions into the discrete 
approximation to the Leray projection. Although not reported here, we 
found that in various experiments this led to greater accuracy than simply 
using normal boundary conditions. 

• The RBF-based Leray projection method only requires boundary condi¬ 
tions on the velocity and does not require any pressure boundary condi¬ 
tions. 

• The pressure can be obtained directly from the discrete velocity using 
(27); solving a Poisson equation for the pressure is unnecessary. 
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• Temporal accuracy is governed purely by the time integrator as the method 


does not split the equations in time. 

• High order spatial accuracy is possible for both the velocity and the pres¬ 
sure. 

The method currently incurs two dimension-dependent costs: an initialization 
cost of and a cost of 0{(PN^) per time-step, where N is the total 

number of nodes. This is a consequence of the use of global RBFs to approxi¬ 
mate the Leray projector. Ongoing work includes reducing the computational 
cost of constructing the discrete Leray projector to 0{dN) using a partition of 
unity type approach and combining this with fast RBF-based methods for ap¬ 
proximating spatial derivatives such as RBF-FD [12, Ch. 5] or RBF-PUM [32]. 
Our hope is that this will enable efficient 2D and 3D simulations on large node 
sets. Another goal is to extend the method to the solution of the incompressible 
Navier-Stokes equations on irregular 2D and 3D domains. This will entail the 
development of robust strategies for meshfree advection that work in tandem 
with the RBF-based Leray projection. 
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